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SUMMARY 


Collection efficiency and ice accretion calculations were made for a commercial transport 
using the NASA Lewis LEWICE3D ice accretion code, the ICEGRID3D grid code and the 
CMARC panel code. All of the calculations were made on a Windows 95 based personal com- 
puter. The ice accretion calculations were made for the nose, wing, horizontal tail and vertical tail 
surfaces. Ice shapes typifying those of a 30 minute hold were generated. Collection efficiencies 
were also generated for the entire aircraft using the newly developed unstructured collection effi- 
ciency method. The calculations highlight the flexibility and cost effectiveness of the 
LEWICE3D, ICEGRID3D, CMARC combination. 


NOMENCLATURE 


AAOA 

Aircraft-angle-of-attack, degrees 

d 

Droplet diameter, pm 

HTC 

Convective heat transfer coefficient, W/m l 2 /K 

LWC 

Liquid Water Content, g/m 2 

MVD 

Median Volume Diameter, pm 

W 

Mass flow, kg/s 

P 

Collection efficiency 

a 

Angle-of-attack, degrees 

e 

Radial angle around inlet lip measured from the upper inlet lip, degrees 
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I. INTRODUCTION 


The design of ice protection systems for new and derivative aircraft has typically been an 
expensive and cumbersome procedure. The design work has typically involved a hodgepodge of 
methodologies including 2-Dimensional (2D) experimental data, 2D trajectory and ice accretion 
codes, wind tunnel tests and flight tests. Difficulties can arise in using any of the methods. Exper- 
imental data may not exist for the geometry or conditions of interest. Applying 2D analytical 
methods to aircraft wing and tail configurations involves corrections for sweep, taper, up-wash, 
down-wash, twist, etc. many of which are not exact. Wind tunnel tests can involve scaling issues if 
sub-scale models are used or blockage and configuration problems if full scale models are used. 
Where full configuration, full scale effects are important (e.g., fuselage mounted inlets and 
radomes, ducts and flaps) the above methods cannot be employed with any confidence. In these 
cases expensive full scale flight testing has typically been required. 

These problems have led to the development over the last 15 years of a series of 3-Dimen- 
sional (3D) droplet trajectory codes and ice accretions codes (ref. 1-5). These methods involve the 
use of various 3D flow codes and droplet trajectory codes which can aid the designer in determin- 
ing where and to what level ice protection is needed. The increased sophistication of these codes 
along with the decreased cost of computing cycles has allowed the designer to get elegant, cheap 
solutions for complex configurations early on in the design cycle decreasing the overall develop- 
ment costs of these systems. 

This paper outlines what is probably the least expensive methodology available to the 
designer of the ice protection system. The calculations presented herein were made using the 
LEWICE3D ice accretion code (ref. 3), the ICEGRID3D unstructured cartesian grid code (ref. 6), 
and the CMARC 3D panel code (ref. 7-9) on a personal computer. These tools can be used by 
experienced users to generate results for full scale configurations in several days on a Personal 
Computer (PC). 

The configuration and conditions for the study were chosen to exemplify those of a typical 
ice protection system design. Calculations were made for the wing, nose and tail of a commercial 
transport in a hold condition. Results for pressure coefficient, collection efficiency, heat transfer, 
and ice shapes on these surfaces are presented. In addition results for the collection efficiency of 
the entire aircraft using the new unstructured collection efficiency method are presented. 

II. ANALYTICAL METHOD 

The icing calculation required a three step process. The ICEGRID3D code was used to 
generate the grid for the droplet trajectory calculations. The CMARC code was used to generate 
the velocities on this grid and to generate the surface velocities needed in the LEWICE3D code. 
The LEWICE3D code then used the panel model, surface velocity information and the grid to 
make the ice accretion calculations. All of the calculations were made on a a Pentium II PC with a 
single 400 MHZ CPU and 384 MB of RAM. 


2 

American Institute of Aeronautics and Astronautics 



A. ICEGRID3D 


The ICEGRID3D code was developed at Lewis Research Center by Bidwell and Coirier 
(ref. 6) specifically for the task of optimizing droplet trajectory calculations in the LEWICE3D 
code. The program automatically produces grids which are optimal for droplet trajectory calcula- 
tions hence reducing the required effort to produce a “good droplet trajectory grid. The program 
also produces a minimum of grid points which reduces the panel code calculation times. The pro- 
gram requires the surface geometry and an input file describing the grid volume and refinement 
parameters. The code refines the grid near regions of interest which can include, parts of the 
geometry, lines, and points. The code is similar to an oct-tree method (ref. 10) in that it recursively 
divides the original grid volume until the refinement criteria tor each cell has been met. The code 
will not refine cells internal to the geometry. The code is different from most oct-tree methods in 
that the grid volume is allowed to be multiple-skewed, multi-blocked, plus different refinement 
functions can be used in any direction. This last feature is where the code differs from the oct-tree 
methods in that it allows a given cell to be divided into 8, 4 or 2 cells depending on the refinement 
function instead of the oct-tree method which divides a cell into 8 cells if the refinement is 
required. This results in grids with much fewer cells for cases where grid requirements are dispar- 
ate in the different directions (e.g. swept wings which have a much smaller cell size requirement 
in the chord-wise direction than in the spanwise direction). 


B. CMARC 

CM ARC is a first order 3D potential flow panel code (ref. 7) based on the PM ARC (ref. 8) 
and VS AERO codes (ref. 9). Geometries are represented as quadrilaterals which have constant 
doublet and source distribution. The formulation used in CMARC results in a solution that is sec- 
ond order accurate allowing for accurate flow solutions with fewer panels and less CPU time than 
other first order methods. The code can generate solutions for internal and external compressible 
flows and can handle very large problems (-10,000 panels). 


C. LEWICE3D 

The LEWICE3D grid based code incorporates droplet trajectory, heat transfer and ice 
shape calculation into a single computer program (ref. 3). The code can handle generic multi- 
block structured grid based flow solutions, unstructured grid based flow solutions, simple carte- 
sian grids with surface patches, and adaptive grids with surface patches. The latter two methods 
allow the use of generic panel code input which is a computationally efficient method for generat- 
ing ice shapes.The code can handle overlapping and internal grids and can handle multiple planes 
of symmetry. Calculations of arbitrary streamlines and trajectories are possible. The code has the 
capability to calculate tangent trajectories and impingement efficiencies for single droplets or 
droplet distributions. Ice accretions can be calculated at arbitrary regions of interest in either a 
surface normal or tangent droplet trajectory direction. 

The methodology used in the LEWICE3D analysis can be broken into six basic steps for 
each section of interest at each time step. In the first step, the flow field is generated by the user. 


3 

American Institute of Aeronautics and Astronautics 



Second, surface streamlines are calculated. The surface streamline analysis uses a variable step 
size fourth-order Runge-Kutta integration scheme developed by Bidwell (ref. 3). Third, tangent 
trajectories are calculated at the region of interest. An array of particles is released between the 
tangent trajectories in the fourth step. These impacting particles are used to calculate collection 
efficiency as a function of surface position. The droplet trajectory analysis is basically that of Hill- 
yer Norment (ref. 1) with modifications by Bidwell (ref. 3). At the heart of the droplet trajectory 
analysis is the variable step predictor-corrector integration scheme by Krogh (ref. 1 1 ). The fifth 
step involves interpolating or extrapolating the collection efficiencies onto the streamlines. In the 
sixth step the ice accretion for the streamline is calculated. The ice accretion model is basically 
that of the LEWICE2D code applied along surface streamlines (ref. 12,13). 

An unstructured collection efficiency calculation method was added to the LEWICE3D 
code recently and results using the method are presented. This method, which was developed by 
Bidwell, was developed as a simple, efficient means of generating collection efficiency maps for 
large complex surfaces. The user specifies an upstream release region (a simple vertical quadrilat- 
eral) which encompasses the region of interest and a minimum cell size. The region is then recur- 
sively subdivided (trajectories are released at the node points of the subdivided region) in areas 
where particle trajectories become close or hit the body until the minimum cell size requirement is 
met. The collection efficiency for each panel is then equal to the number of particles that hit the 
panel times the area of the minimum upstream cell size divided by the area of the panel. The 
method efficiently incorporates the tangent droplet trajectory search and the collection efficiency 
calculation into one step. Although the method is more computationally expensive than traditional 
methods (many more droplet trajectory integrations are required), it is not prohibitive. Large prob- 
lems can still be handled with overnight runs on a PC (4-20 hours). The unstructured method has 
the advantage that large, complex problems are easier to setup than for more traditional methods. 
Is is also advantageous in that the method can handle complex geometries (inlet lip-duct prob- 
lems, highly separated flows, multi-element) that other more traditional codes may not handle. It 
can also find impingement regions (multi-element) that other codes may miss if not expected. 

III. CONFIGURATION 

The conditions and geometry for the analysis were chosen to typify a commercial trans- 
port in a hold condition. The configuration used in the analysis, which is similar to a Boeing 737 
airplane, was chosen to be consistent with those of a wind tunnel model used at Langley Research 
Center for iced stability and control measurements (ref. 14,15). The coordinates for this model 
were obtained from Langley Research Center. A hold condition was chosen because it was an 
important design point for the ice protection systems on a commercial transport. 

The aircraft model was a simplified version of the Boeing 737 airplane. The engine 
nacelles and engine pylons were not included in the flow model. The wings were modeled using a 
single element model with retracted flaps and slats. The deletion of these features from the aircraft 
should not be deemed too severe, considering the position of the engines and the position of the 
flaps and slats in a typical hold condition. 

A single flight condition was chosen for the study. The flight condition for the analysis 
was; aircraft-angle-of-attack (AAOA), 4 degrees, altitude, 4572 meters, airspeed, 135 m/s, static 
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temperature, 259 K, icing time, 30 minutes, liquid water content (LWC), 0.51 g/m , and median 
volume diameter (MVD), 20 pm. 

The computational parameters were chosen from experience, correlations and the desire to 
limit the computational resources required. A single drop size and icing time step were chosen for 
the study. A LEWICE roughness parameter of 0.5 mm was used for all of the cases (ref. 12). The 
calculations were made at six spanwise stations on the wing, three spanwise stations on the hori- 
zontal and vertical tails, and a single vertical cut on the nose of the aircraft. In the wing analysis, 
three spanwise stations were distributed on the outer portion of the wing and three were distrib- 
uted on the inner portion of the wing. The break in leading edge sweep angle near the 30% span 
location denotes the boundary between the outer and inner portions of the wing. The spanwise 
cuts were located at the 10%, 50% and 90% locations for the inner and outer parts of the main 
wings, the horizontal tail, and the vertical tails. The 10% and 90% span stations were chosen as 
conservative limits for the LEWICE3D methodology. Outside of these limits the spanwise gradi- 
ents can be large and the strip assumptions built into the LEWICE3D methodology can be vio- 
lated. 


IV. ANALYSIS 

Surface velocity, heat transfer, collection efficiency, and ice shapes results are presented 
for the commercial airplane at various stations around the wing, horizontal tail, vertical tail and 
fuselage nose. Ice shape calculations were made for a singe condition representing a hold condi- 
tion. Discussions of the panel model, grids, the LEWICE3D program parameters, and of the indi- 
vidual analysis are given below. 

Figure 1 shows the CMARC panel model and the surface cuts for the analysis. A half- 
plane model was employed because a zero degree yaw angle was used in the analysis. The panel 
model contained 6545 panels and required 54 minutes of CPU time on the 400 MHZ Pentium II 
PC to solve. These panels were concentrated near the leading edge of the wing, horizontal and 
vertical tails and the aircraft nose to generate a “good” ice shape calculation. To obtain a reason- 
able representation of the ice shapes approximately 15 to 20 chord-wise panels are required 
within the impingement limits. This constraint resulted in panel models for the wing, and horizon- 
tal and vertical tails which had approximately 100 chord-wise panels. 

Non-standard settings for the CMARC parameters RFF, RCORES, and RCOREW were 
used to generate the off-body velocities. Values of 25., 0.0002 and 0.0002 were used respectively 
(ref. 7). A relatively large value of RFF was used because of the very small panels used near the 
leading edge of the wing, the tail, and the nose of the aircraft. Very small values of RCORES and 
RCOREW (essentially zero) were used to bypass the singularity “fix” used in the velocity subrou- 
tines because it did not work very well. The singularities were handled in a more consistent man- 
ner in the post processing program VELCOND developed by Bidwell. This program replaces the 
off-body velocities near the surface and internal to the surface with surface velocities interpolated 
from the panel center points. This velocity “conditioning” was key in generating smooth and con- 
sistent results. 


5 

American Institute of Aeronautics and Astronautics 



The off-body velocities were generated using the CMVSCAN post processor developed 
by AeroLogic (ref. 7). The post processor used the flow-field solution generated by CMARC to 
generate velocities at the ICEGRID3D grid points. Separate calculations were done for each of 
the ICEGRID3D blocks. These results were combined into a single multi-block grid file for input 
into LEWICE3D. The off-body calculations for the 3,107,543 points required approximately 45 
hours of CPU time on the 400 MHZ Pentium II PC. 

The grid blocking structure is shown in fig. 2a-c. A total of 9 ICEGRID3D grid blocks 
were used in the grid. Each of the blocks were calculated separately and combined before input 
into the LEWICE3D program. These 9 blocks required approximate 48 hours of CPU time on the 
400 MHZ Pentium II PC to generate. The multi-block structure was used to minimize the total 
number of points used, the time required to re-work a grid if necessary, and the overall 
ICEGRID3D calculation time. The multiple-skewed feature of the ICEGRID3D grid blocking 
scheme allows the grid block to be aligned with the dihedral and sweep of the wing or tail to min- 
imize the number of grid points required, but for cases with large sweep or dihedral change or ele- 
ments with different sweep a single grid block is not optimal. In these cases, multiple blocks were 
used with separate blocks tailored for each element of interest. The use of separate blocks also 
reduces the time to generate a new grid if a region needs to be re-worked because only the grid 
block associated with the region of interest needs to be recalculated. Breaking the problem into 
smaller blocks also reduces the memory required for ICEGRID3D and allows the use of multiple 
machines (i.e. one for each grid) to decrease the turnaround time. 

Figure 2 reveals the overall blocking strategy. One can see that the grid blocks are overlap- 
ping and have various sizes, and skewness. Three grid blocks were used for the wing to handle the 
changes in sweep and dihedral of the wing and to handle the transition into the fuselage grid 
blocks. Transition grids between airfoils and the fuselage were needed for the wing and tails sec- 
tions because the airfoil surfaces had different grid resolution requirements then the fuselage. The 
wings and tails required very few span-wise planes but required many chord-wise points and 
points normal to the surface. The fuselage required the same grid spacing in all three directions 
and hence typically required more points. The three wing grid blocks enveloped the area in front 
of the wing to just aft of the wing. A single grid block was used for the outboard portion of the 
horizontal tail. The vertical tail required two blocks to accommodate the change in sweep. A sin- 
gle grid block was used for the fuselage grid and the transitions from the horizontal and vertical 
tail grids. This fuselage grid block and all of the tail blocks extended from aft of the wing to aft of 
the plane. Two blocks were used to handle the remainder of the fuselage; a block extending from 
the aircraft wind screen to forward of the aircraft, and a block extending from aft of the wing for- 
ward to the wind screen. 

The grid resolution for various regions around the aircraft are shown in fig. 3. The maxi- 
mum cell size for the grids was 0.813 m in all three coordinate directions. For the nose of the air- 
craft a minimum cell size of 0.003 m was used. The same cell size requirements were used in the 
all three coordinate directions for the nose. For the horizontal and vertical tail grids the minimum 
chord-wise and surface normal spacing were 0.003 m. The minimum spanwise spacing for the 
horizontal and vertical tail was 0. 102 m. For the wing the minimum spanwise spacing was 0.203 
m. The wing had a minimum chord-wise and surface normal spacing of 0.002 m. A minimum 
cells size of 0.051 m was used in all three coordinate directions for the mid-fuselage block. 
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The philosophy behind the surface paneling and grid generation is to pack panels and 
points where the trajectories are going to hit the surface and to use few panels and points away 
from the impact regions. The grid cell sizes should be consistent with the panel sizes in the impact 
regions. It does no good to put a lot of grid points in a region with poor paneling or vice-versa. In 
cases where the particles travel close to the body for extended lengths (e.g. the underside of the 
wing at higher angles-of-attack) or where trajectories need to be calculated accurately along the 
fuselage to some region of interest (e.g. radome) smaller panels and grid cell sizes may be 
required. 

Figure 3 shows this philosophy applied to the different elements of the aircraft. For the 
wing, impingement was expected on the underside of the wing due to the configuration of the 
AAOA. Grid points and panels were packed around the leading edge and underside of the wing. A 
small grid spacing was used for a large portion of underside of the wing to carry the tangent tra- 
jectories aft accurately without hitting the body erroneously. For the vertical and horizontal tail, 
the impingement regions were more towards the leading edge and more symmetric due to the 
lower effective angle-of-attack; hence the grids points, and surface panels are packed more sym- 
metrically around these regions. A relatively large spacing was used for the mid-fuselage grid 
because no impingement was expected. The nose grid reflects the impingement expected on the 
nose of the aircraft. The blunt nature of the nose tends to throw particles out away from the body 
if they do not hit, resulting in a shadow zone around the entire fuselage aft of the nose. For this 
reason only a small region around the nose was refined. 


The icing analysis for the commercial aircraft are shown in fig. 4-23. The section profiles 
and coefficient of pressure distributions were generated along the flow directions. The collection 
efficiency, heat transfer coefficient and ice shape profiles were generated in a direction normal to 
the leading edge. Each section-of-interest required approximately 2 minutes of CPU time on the 
400 MHZ Pentium II PC. The overall collection efficiency results shown in fig. 23 involved the 
calculation of approximately 340,000 trajectories and required about 25 hours of CPU time. 

The computational results for the wing are shown in fig. 4-8. Figure 4 depicts the airfoil 
section at the six spanwise stations where the analysis was done. In general the wing was com- 
plex, having taper, twist and section variation. 

Figure 5 shows the coefficient of pressure distribution at each of the six spanwise stations. 
The section lift coefficient increases as the wing root is approached on the outboard section of the 
wing due to wing twist and up-wash effects. The effective angle-of-attack for the inner five sec- 
tions is relatively constant. The leading edge coefficient of pressure gradient decreases (i.e. the 
peaks are more rounded) as the root is approached due to the increase in chord and sweep. 

Figure 6 shows the heat transfer coefficient distribution results for the wing. In general, as 
the wing root was approached the stagnation heat transfer and the peak heat transfer on both the 
upper and lower surfaces decreased. Additionally the laminar-turbulent transition points on the 
upper and lower surfaces moved further aft as the spanwise location moved towards the wing root. 
These effects are due to the increase in chord and sweep angle as the root was approached and the 
corresponding decrease in leading edge velocity gradients (and the associated coefficient of pres- 


7 

American Institute of Aeronautics and Astronautics 



sure). The heat transfer coefficient distribution is calculated from the surface velocity distribution 
in the integral boundary layer technique employed, hence, trends seen in the heat transfer distribu- 
tion can be gleaned from the coefficient of pressure distribution (i.e where there are large gradi- 
ents in the coefficients of pressure there are high heat transfer coefficients). 

The collection efficiency results for the wing are shown in fig. 7. In general the shapes of 
the collection efficiency curves were similar except for the two most inboard sections. The shapes 
are similar for the four outboard sections because the airfoil sections employed at these locations 
were similar. The inner two sections have blunter sections, hence the collection efficiency curves 
have wider more rounded peaks. The curves show that the majority of the water is hitting on the 
underside of the wing. This was due to the high AAOA (4 degrees) used in the study. In general 
the maximum collection efficiency decreased as the root was approached because of the increase 
in chord of the wing. 

The calculated ice shapes for the wing are shown in fig. 8. Due to the mixed or glaze 
nature of the ice (i.e. there was significant run-back at the leading edge) the ice shapes follow the 
trends set by the heat transfer coefficient described above. When there is more water impinging 
than can freeze at the leading edge, the amount of ice formed is determined by the heat transfer at 
the leading edge. The result being that the ice shape thickness distribution and the resulting ice 
shape take on characteristics of the heat transfer coefficient distribution. From the figures it can be 
seen that the ice thickness at the stagnation zone decreases as the root is approached due to the 
decreasing stagnation heat transfer. Also the maximum thickness on the upper and lower surfaces 
moves further aft as the root is approached due to the aft movement of the boundary layer transi- 
tion point. 

Results for the horizontal tail are presented in fig. 9-13. The airfoil sections for the tail are 
shown in figure 9. The horizontal tail element was of a more simpler construction than the wing, 
having no twist and constant sweep, section and taper ratio. 

Figures 10, 1 1 show the coefficient of pressure and heat transfer coefficient distributions 
respectively. From the coefficient of pressure distribution a decrease in angle -of-attack can be 
seen as the root is approached due to down-wash. This results in a slight shift in the stagnation 
region, and hence heat transfer coefficient distribution, towards the upper surface. Also evident is 
a decrease in the leading edge pressure gradient due to the increase in chord. This results in a 
decrease in the stagnation heat transfer and the maximum heat transfer on both the upper lower 
surfaces as the root is approached. 

Collection efficiency distributions are shown in fig. 12 for the horizontal tail.The curves 
have similar shapes and show a decrease in maximum collection efficiency as the wing root is 
approached due to the increase in chord. Also evident is a shift in the impingement region towards 
the upper side surface of the horizontal tail. This can be attributed to the decrease in angle-of- 
attack as the root is approached, due to the wing down-wash. 

The calculated ice shapes for the horizontal tail are shown in fig. 13. The results for the 
horizontal tail are similar to those for the wing. From fig. 13, a decrease in the stagnation region 
ice thickness from tip to root can be seen due to the decrease in stagnation heat transfer. Also, the 
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maximum thickness on the upper and lower surfaces moves further aft as the root is approached 
due to the aft movement of the boundary layer transition point. There is also a slight shift in the 
ice shape towards the upper surface as the root is approached due to the decreased angle-of-attack 
associated with the wing down-wash. 

Figures 14-18 depict the icing analysis for the vertical tail. The vertical tail airfoil sections 
are shown in fig. 14. The vertical tail was of a simple construction containing symmetric sections 
of a similar airfoil type which varied in chord length due to the taper and the change in leading 
edge sweep angle. 

The coefficient of pressure and heat transfer coefficient distributions for the vertical tail 
are shown in fig. 15 and 16. As for the wing and the horizontal tail the leading edge pressure coef- 
ficient of pressure gradient decreases and the stagnation heat transfer coefficient and maximum 
heat transfer coefficient on the upper and lower surface decrease as the as the fuselage is 
approached. This is due to the increase in chord and sweep angle as the fuselage is approached. 
Also noteworthy is a slight discontinuity in the pressure distribution near the leading edge of the 
outer two spanwise stations. This is due to the coarse geometry definition of the original model. 
Many points were splined into the leading edge region from a very few original leading edge 
points (3-5 points) and spanwise points. This can cause wrinkles and slope discontinuities which 
show up as discontinuities in the pressure distribution. 

Collection efficiency distributions for the vertical tail are shown in fig. 17. The maximum 
collection efficiency decreases and the impingement limits increase as the fuselage is approached 
due to both the increase in chord and the increase in sweep angle. The maximum collection for the 
outer two sections-of-interest were closer than expected. This was due to a deficiency in the grid 
density at the outer portion of the vertical tail. Evidence of insufficient grid spacing include flat- 
tened or in extreme cases dips in the collection efficiency curve near the maximum. 

The ice shapes for the vertical tail are shown in fig. 18. The resulting ice shapes show the 
same trends as the wing and horizontal tail. The maximum ice thickness decreases and the point at 
which the maximum ice thickness occurs moves further aft as the wing root is approached. 

The icing analysis for the fuselage nose is shown in fig. 19-22. Only a single vertical ice 
shape was calculated for the nose. Other radial cuts at the nose could have been generated but the 
results would have been similar. The nose of the aircraft was very blunt, much like a large sphere. 
This led to a very small leading edge coefficient of pressure gradient and hence to small values of 
heat transfer coefficients in the stagnation region. The maximum collection efficiency for the nose 
was also very small (-.20) due to the large, blunt nature of the nose. The ice shape which is shown 
in fig. 22 is very thin with the maximum thickness point being relatively far back on both the 
upper and lower surfaces. This is due to the very low stagnation heat transfer and the large laminar 
flow region at the leading edge of the nose. 

Figures 23a-f depict the overall picture of the collection efficiency for the airplane. The 
collection efficiency for the aircraft was generated using the new unstructured collection effi- 
ciency method developed by Bidwell. The calculations, although computationally expensive, give 
a much broader picture of the collection efficiency characteristics of the airplane. Fig. 23a shows 
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the aircraft along with the unstructured droplet trajectory release point grid used to generate the 
collection efficiency map. The grid shows the release points of the trajectories used in the unstruc- 
tured calculation, hence, it is located approximately an aircraft length ahead of the body. Each of 
the grid nodes represent a particle release point used in the unstructured droplet trajectory itera- 
tion. The methodology involves recursively refining the grid around trajectories that are close to 
or hit the body until a minimum cell size is attained. The result being that the grid is very dense 
around regions that correspond to areas where the trajectories hit the body. The grid was superim- 
posed on the image to give the analyst insight into the shape and size of the free-stream impinge- 
ment region. 

The free-stream wing impaction region is the long region which extends from the center- 
line outboard (fig. 23a). The changes in effective angle-of-attack and wing dihedral can be seen in 
the bowed shape of the wing impingement region. For a fixed section, higher angles-of-attack will 
result in a downward movement of the impacting droplet trajectory release points. For a fixed 
angle-of-attack, a higher section will result in an upward movement of the impacting droplet tra- 
jectory release points. From fig. 23a, the maximum effective angle-of-attack of the wing can be 
seen to occur somewhere around the 30% span location. From fig. 23c the increase in maximum 
collection efficiency can be seen as the root is approached due to the increase in chord. Also note- 
worthy is the drop-off in collection efficiency at the wing root due to “shadowing” by the nose of 
the airplane. 

The free-stream impaction region for the horizontal tail is located above the wing impac- 
tion region (fig. 23a). It is shorter and bowed downward. The region shows the spanwise changes 
in dihedral and effective angle-of-attack. From fig. 23a it can be seen that the maximum angle-of- 
attack for the horizontal tail occurs near the tip. This trend was consistent with that observed for 
the coefficient of pressure. Similar to the wing case, it can be seen that the maximum collection 
efficiency drops off as the root is approached due to the increase in chord and that there is a 
“shadow” zone near the root caused by fuselage (fig. 23d). 

The results for the vertical tail are shown in figures 23a and 23e. The free-stream impac- 
tion region for the vertical tail is vertical and occurs along the symmetry plane because the aircraft 
flow was calculated for a zero yaw angle. From fig. 23e, the fall-off in collection efficiency as the 
root is approached due to increasing chord and the absence of impingement near the root of the 
vertical tail due to fuselage “shadowing” can be seen. 

Figures 23a and 23f depict the results for the aircraft nose. The free-stream impaction 
region and the impingement region are almost circular due to the blunt spherical nature of the air- 
craft nose. Because of the size and bluntness of the nose the collection efficiencies are much lower 
than for the wing and tail elements. 


V. CONCLUSION 

The grid based LEWICE3D, CMARC, ICEGRID3D combination proved to be an inex- 
pensive, flexible, accurate ice protection system design tool. The flow and ice accretion calcula- 
tions were done quickly, and inexpensively for a range of aircraft elements. 


10 

American Institute of Aeronautics and Astronautics 



The ICEGRID3D, CMARC, grid based LEWICE3D combination was, in general, inex- 
pensive to operate for the full aircraft configuration. The grid calculation took approximately 45 
hours and the CMARC flow calculations took about 50 hours on the 400 MHZ Pentium II PC 
computer. The ice accretion calculations required about 2 minutes per section-of-interest. The 
overall collection efficiency calculations took approximately 25 hours. 

In general, the results for coefficient of pressure, heat transfer, collection efficiency and ice 
shape appear reasonable. The coefficient of pressure followed the correct trends for sweep. The 
stagnation heat transfer coefficient and the maximum heat transfer coefficient on the upper and 
lower surfaces decreased with decreasing leading radius of curvature as expected. The location of 
the maximum heat transfer coefficient also moved further aft with decreasing leading edge radius. 
The maximum collection efficiency decreased with increased chord and increased sweep angle. 
Because of the mixed nature of the icing condition used in the study, the ice shapes followed 
trends set by the heat transfer coefficient. The stagnation ice thickness decreased with increased 
chord (i.e. decreased leading edge radius). The position of the maximum thickness for both the 
upper and lower surfaces moved further aft with decreased leading edge radius. 

The overall collection efficiency plots generated using the new unstructured collection 
efficiency method, although relatively computationally expensive, gave greater insight into the 
icing characteristics of the airplane. From the overall plots the effects of sweep, effective angle- 
of-attack and chord can be seen more easily than using traditional cuts. The “shadowing” effects 
of the fuselage on the wings and tail are also easily seen from the overall collection efficiency 
maps. 


II 

American Institute of Aeronautics and Astronautics 



REFERENCES 


1 . Norment, H.G., “Calculation of Water Drop Trajectories To and About Three-Dimensional 
Lifting and Non-lifting Bodies In Potential Airflow,” NASA CR-3935, October 1985. 

2. Kim, J J., “Particle Trajectory Computation on a 3-Dimensional Engine Inlet ” NASA CR- 
175023, January 1986. 

3. Bidwell, C.S., and Potapczuk, M.G., “Users Manual for the NASA Lewis Three-Dimen- 
sional Ice Accretion Code (LEWICE3D),” NASA TM- 105974, December 1993. 

4. Nathman, J.K., ”ICE Particle Trajectory and Ice Accretion Program, User’s Manual, Ver- 
sion 2.0,” December 1992. 

5. Farrel, D.M., Vittal, B. V.R., “Computer Simulation of Water Ingestion for the T800-LHT- 
801 Engine Anti-Icing Analysis,” AIAA-96-2652. 

6. Bidwell, C.S., “Collection Efficiency and Ice Accretion Calculations for a Boeing 737-300 
Inlet,” NASA TM- 1 07347, October 1 996. 

7. Pinella, D., Garrison, P., “CM ARC, Three-Dimensional Low Order Panel Code, Version 
2.1 ” 1997. 

8. Ashby, D.L., Dudley, M.R., Iguchi, S.K., “Development and Validation of an Advanced 
Low-Order Panel Method,” NASA TM- 101024, October 1988. 

9. Maskew, B., “Program VSAERO Theory Document: A Computer Program for Calculating 
Nonlinear Aerodynamic Characteristics of Arbitrary Configurations, ” NASA CR-4023, 
September 1987. 

10. Johnson, F.T., Samant, S.S., Bieterman, M.B., Melvin, R.G., Young, D.P., Bussoletti, J.E., 
and Hilmes, C.L., “TranAir: A Full-Potential, Solution-Adaptive, Rectangular Grid Code 
for Predicting Subsonic, Transonic, and Supersonic Flows About Arbitrary Configura- 
tions,” NASA CR-4348, December 1992. 

1 1 . Krogh, F.T., “Variable Order Integrators for Numerical Solutions of Ordinary Differential 
Equations,” Jet Propulsion Lab Technology Utilization Document No. CP-32308, Novem- 
ber 1970. 

12. Ruff, G.A., Berkowitz, B.M., “Users manual for the NASA Lewis Ice Accretion Prediction 
Code (LEWICE) ,” NASA CR- 185 129, May 1990. 

13. Wright, W.B., “Users manual for the Improved NASA Lewis Ice Accretion Code LEWICE 
1.6,” NASA CR- 198355, Junel995. 


12 

American Institute of Aeronautics and Astronautics 



14. Reehorst, A., Potapczuk, ML, Ratvasky, T., and Gile Laflin, B., “Wind Tunnel Measured 
Effects on a Twin-Engine Short-Haul Transport Caused by Simulated Ice Accretions,’ 
AIAA-96-0871, NASA TM- 107 143, January 1996. 

15. Reehorst, A., Potapczuk, M., Ratvasky, T., and Gile Laflin, B., “Wind Tunnel Measured 
Effects on a Twin-Engine Short-Haul Transport Caused by Simulated Ice Accretions, Data 
Report,” NASA TM- 1074 19, May 1997. 


13 

American Institute of Aeronautics and Astronautics 



Figure 1 . - CMARC panel model and sections-of-interest. 



(a) Front view. 


Figure 2. - CMARC panel model and 1CEGRID3D grid block structure. 


14 

American Institute of Aeronautics and Astronautics 


(b) Side view. 



(c) Bottom view. 

Figure 2. - Concluded. CMARC panel model and ICEGRID3D grid block structure. 


15 

American Institute of Aeronautics and Astronautics 




(e) Side view of fuselage. 

Figure 3. - CMARC panel model and ICEGRID3D grid structure. 
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Figure 4. - Wing profile at several span-wise stations. 



AXIAL DISTANCE (NONDIMENSIONAL) 


Figure 5. - Coefficient of pressure for wing at several span-wise stations for aircraft-angle-of- 
attack, 4 degrees. 
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Figure 6. - Analytical heat transfer coefficient distribution for the wing at several span-wise sta- 
tions. Flight conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 
20 pm; static temperature, -14° C; static pressure, 57800 Pa. 



Figure 7. - Analytical collection efficiency distribution for the wing at several span-wise stations. 
Flight conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 pm; 
static temperature, -14° C; static pressure, 57800 Pa. 
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Figure 8. - Calculated ice accretions for the wing at several span-wise stations. Flight conditions; 
airspeed, 1 35 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 pm; static tempera- 
ture, -14° C; static pressure, 57800 Pa. 



Figure 9. - Horizontal tail section profiles at several span-wise stations. 
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Figure 10. - Coefficient of pressure for horizontal tail at several span-wise stations for aircraft- 
angle-of-attack, 4 degrees. 
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Figure 11.- Analytical heat transfer coefficient distribution for the horizontal tail at several span- 
wise stations. Flight conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; 
drop size, 20 pm; static temperature, -14° C; static pressure, 57800 Pa. 
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Figure 12. - Analytical collection efficiency distribution for the horizontal tail at several span-wise 
stations. Flight conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop 
size, 20 pm; static temperature, -14° C; static pressure, 57800 Pa. 



Figure 13. - Calculated ice accretions for the horizontal tail at several span-wise stations. Flight 
conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 pm; static 
temperature, -14° C; static pressure, 57800 Pa. 
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Figure 14. - Vertical tail section profiles at several span-wise stations. 



0.00 0.10 0.20 0.30 0 40 0.50 0.60 0.70 0.80 0.90 1.00 


AXIAL DISTANCE (NONDIMENSIONAL) 

Figure 15. - Coefficient of pressure for vertical tail at several span-wise stations for aircraft-angle- 
of-attack, 4 degrees. 
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Figure 16. - Analytical heat transfer coefficient distribution for the vertical tail at several span- 
wise stations. Flight conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees, 
drop size, 20 |im; static temperature, -14° C; static pressure, 57800 Pa. 
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Figure 17. - Analytical collection efficiency distribution for the vertical tail at several span-wise 
stations. Flight conditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop 
size, 20 pm; static temperature, -14° C; static pressure, 57800 Pa. 
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Figure 18. - Calculated ice accretions for the vertical tail at several span-wise stations. Flight con- 
ditions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 pm; static 
temperature, -14° C; static pressure, 57800 Pa. 
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Figure 19. - Coefficient of pressure for fuselage nose for aircraft-angle-of-attack, 4 degrees. 
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Figure 20. - Analytical heat transfer coefficient distribution for the fuselage nose. Flight condi- 
tions; airspeed, 1 35 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 pm; static tem- 
perature, -14° C; static pressure, 57800 Pa. 
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Figure 21. - Analytical collection efficiency distribution for the fuselage nose. Flight conditions; 
airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 pm; static tempera- 
ture, -14° C; static pressure, 57800 Pa. 
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Figure 22. - Calculated ice accretions for the fuselage nose. Flight conditions; airspeed, 135 m/s; 
aircraft-angle-of-attack, 4 degrees; drop size, 20 p.m; static temperature, -14° C; static 
pressure, 57800 Pa. 



Figure 23. - Unstructured collection efficiency calculation for aircraft. Flight conditions; airspeed, 
135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 Jim; static temperature, -14° C; 
static pressure, 57800 Pa. 
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Figure 23. - Concluded. Unstructured collection efficiency calculation for aircraft. Flight condi- 
tions; airspeed, 135 m/s; aircraft-angle-of-attack, 4 degrees; drop size, 20 p.m; static tem- 
perature, -14° C; static pressure, 57800 Pa. 
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